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ABSTRACT 

Most previous calculations of the annihilation decay rates of heavy quarkonium 
were based on factorization assumptions that were unproven and, in some cases, 
even incorrect. The recent development of a general factorization formula for heavy 
quarkonium annihilation rates has provided a rigorous theoretical foundation for 
such calculations. The factorization formula is based on the use of the effective field 
theory NRQCD to factor the decay rate into short-distance coefficients that can be 
calculated in perturbation theory and long-distance matrix elements that can be 
computed using lattice simulations. This approach allows annihilation decay rates 
to be computed entirely from first principles. 



1. History of Quarkonium Annihilation Calculations 

1.1. Leading-order Calculations (1974-1976) 

The discovery of charmonium in November 1974 launched a revolution in particle 
physics. Within days of the announcement of the discovery, there were papers 
submitted to Physical Review Letterstl which interpreted the new resonances as 
bound states of a charm quark and antiquark analogous to positronium. Among 
other things, the authors calculated the annihilation decay rates of the S-wave 
states of charmonium, which were eventually named rjc and J/ip. They assumed 
that charmonium decays into light hadrons via the annihilation of the cc pair into 
gluons, just as positronium decays by annihilation into photons. The decay rates for 
positronium are proportional to |i?5(0)p, where Rs{0) is the radial wavefunction at 
the origin. The decay rates for rjc ^ gg, rjc ^ •yy, ip ^ ggg, and tfj e~^e~ can all 
be obtained to leading order in as from the corresponding positronium calculations 
by replacing a by as and rUe by M (where M is the charm quark mass) in the 
coefficient of |i?5(0)p and multiplying by an overall color factor. 

The annihilation decay rates of P-wave charmonium states were calculated to 
leading order in ag by Barbieri and coIIaboratorg^'H in 1976. The decay rates were 
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assumed to be proportional to |i?p(0)p, where -Rp(O) is the derivative of the radial 
wavefunction at the origin. This assumption seemed to be valid for the decay rates 
of the spin-0 and spin-2 states into light hadrons via Xco ~^ 99 and Xc2 99- It 
also seemed to be valid for their electromagnetic decays Xco 11 and Xc2 11- 
However, for decays of the spin-1 states into light hadrons, the partial decay rates 
for he — * 999 and Xa <l<i9 were infrared divergent. These divergences were 
interpreted as a signal that the decay rates are sensitive to the binding energy and 
therefore to nonperturbative effects of QCD. 

In 1976, the ITEP groupQ wrote a review of applications of QCD to charmonium, 
including the calculations of annihilation decay rates. While there was dramatic 
improvement in experimental measurements of the decay rates in subsequent years, 
their theoretical discussion of charmonium decay rates remained for the most part 
up-to-date for over 15 years. Implicit in the calculations of the decay rates is the 
assumption that they can be factored into a nonperturbative part involving the 
wavefunction and a perturbative part, which can be calculated in terms of the 
running coupling constant as{M). 

1.2. Next-to-leading- order Calculations (1979-1981) 

The first calculation of the radiative corrections for quarkonium decays was 
carried out in 1979 by Barbieri et alB. They calculated the next-to-leading order 
corrections to the decay rates for 77^ 77 and rj^. LH, where LH represents light 
hadrons. The next-to-leading order radiative correction for ip LH was calculated 
by Mackenzie and Lepage§ in 1981. These calculations were based on the assump- 
tion that the decay rate could be factored into [/^^(O)^ multiplied by a perturbative 
coefficient. There was no general theory guaranteeing such a factored form, but the 
results of these calculations provide empirical evidence that the assumed factoriza- 
tion occurs, at least to next-to-leading order in as- 

The next-to-leading order corrections to the decay rates of the P-wave states 
Xco and Xc2 were calculated by Barbieri and collaborators! in 1980. The decay 
rates were assumed to factor into |i?p(0)p multiplied by a perturbative coefficient. 
For the electromagnetic decays Xco 11 and Xc2 11, the next-to-leading order 
corrections were consistent with the factorization assumption. However, for the 
inclusive decays into light hadrons Xco LH and Xc2 — ^ LH, the next-to-leading 
order correction was infrared divergent. As in the case of the decays of he and xa 
at leading order, these divergences were interpreted as a signal that the decay rates 
were sensitive to nonperturbative effects. 

1.3. The Breakthrough (1 992- 1 994 ) 

In the decade following the calculations of the radiative corrections for the S- 
wave and P-wave states, there was very little progress in the theory of heavy quarko- 
nium annihilation. The breakthrough came in 1992 with the solution of the problem 
of infrared divergences in the P-wave decay ratesH. The breakthrough was stimu- 
lated in part by Experiment E760 at Fermilab, which produced charmonium states 
via resonant pp-annihilation and measured their masses, widths, and branching 



fractions with unprecedented precisionE The breakthrough was made possible by 
the development of nonrelativistic QCDt3, an effective field theory for bound states 
of heavy quarks and antiquarks. The NRQCD approach has taken us far beyond 
the solution to the P-wave problem. It provides a general factorization formula for 
annihilation decay rates of heavy quarkonium, including not only perturbative cor- 
rections but relativistic corrections as welltll. It places the entire subject of heavy 
quarkonium annihilation on a firm theoretical foundation. 

2. Problems in the Conventional Approach 

2.1. S-wave Decays 

Previous calculations of the annihilation decay rates of heavy quarkonium were 
based on the assumption that the decay rate factors into a short-distance part, 
which can be calculated as a perturbation series in the running coupling constant 
as{M), and a long-distance factor, which depends on the nonperturbative dynamics 
of the bound state. To illustrate how such calculations are carried out in practice, 
we review the calculation of the decay rate of the rjc to next-to-leading order in as- 

The 7]c is predominantly a |cc) state, and it decays into light hadrons through 
the annihilation of the cc pair. The annihilation rate depends on the wavefunction 
of the cc pair, which is essentially nonperturbative. Instead of trying to compute 
the decay rate of the bound state directly, we calculate the decay rate of a simpler 
system in perturbation theory. This system consists of a c and c that are above 
threshold and on their mass-shells, with a prescribed nonrelativistic wavefunction. 
The wavefunction is color-singlet (^^■'/v^), spin-singlet ((ti — l'\)/y/2), and S-wave 
{ip{r) = i?(r)/-\/47r, with / r'^\R{r)\'^dr = 1). At leading order in perturbation theory, 
the cc annihilates into 2 gluons. At next-to- leading order in as, we must add the 
rates for cc ggg and cc — > qqg, as well as the interference term between the 
amplitudes of order as and a^ for cc gg. Ultraviolet divergences are removed by 
renormalization of the coupling constant. Taking the nonrelativistic limit, we find 
that the decay rate is proportional to the square of the radial wavefunction at the 
origin i?(0): 
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The term is the average value of the inverse of the velocity when weighted by 

the wavefunction. This term, which comes from Coulomb scattering of the c and 
c, diverges in the nonrelativistic limit. Like -R(O), it is sensitive to nonperturbative 
effects involving distances much longer than 1/M. We can factor the decay rate (|l|) 
into long-distance and short-distance parts as follows: 

r(c.) ^ \Rm^ (i + =|i(iA>) X 2^ (i + 51^^) p) 

We now return to the rj^ which also decays via the annihilation of the cc pair. 
We assume that this decay rate can also be factored into short-distance and long- 



distance parts and that the short-distance part is insensitive to whether the c and 
c are shghtly above threshold and on their mass sheUs, as in the calculation above, 
or slightly below threshold in a bound state. The decay rate is then given by an 
expression analogous to (j^), with the same short-distance factor but a different 
long-distance factor: 
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This next-to-leading order result was first obtained by Barbieri et ala in 1979. 
The correction (7ra5/3)(l/f), which arises from Coulomb scattering of the cc pair, 
has been factored into -R5(0). This term is universal, appearing also in the next- 
to-leading order corrections for rjc 77, ip LH, and ip e^e^ . Thus one can 
consistently equate the long-distance factor |_R5(0)p in all of the S-wave decay rates. 

2.2. P-wave Decays 

Given the success of the factorization assumption for S-wave decays, it is natural 
to assume that P-wave decays factor in a similar way. Instead of |-R5(0)p, the cor- 
responding long-distance factor would be |i?p(0)p, where -Rp(O) is the derivative of 
the wavefunction at the origin. Using this factorization assumption, the decay rates 
of the four P-wave states into light hadrons were calculated by Barbieri and various 
collaboratordiffllli to order a^. In the case of Xco and Xc2, this is next-to- leading 
order in as- The next-to-leading order corrections include a Coulomb singularity 
{4:nas/S){l/v) that can be factored into |i?p(0)p. The resulting expressions for the 
decay rates have the form 
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where Uf = 3 is the number of light quark fiavors, the ... represent numerical 
constants, and A is an infrared cutoff on the energies of gluons in the final state. 
Logarithmic infrared divergences appear as we take the limit A — »• 0. It is clear 
from the expressions in (^-(|^) that the infrared divergences can not be factored 
into |-Rp(0)p. In phenomenological applications, the infrared cutoff A has been 
identified with the inverse of the radius of the bound state, or with its binding 
energy, or with the inverse of the confinement radius. It is evident, however, that 
the presence of infrared divergences implies a failure of the factorization assumption 
upon which the entire calculation is based. 



2.3. Questions 



The explicit calculations of quarkonium annihilation rates carried out from 1974 
to 1981 raise a number of important questions that must be answered by a rigorous 
theory of annihilation. 

1. Do the annihilation rates factorize? Empirically, factorization holds to next- 
to-leading order in for S-waves and also for electromagnetic decays of P-waves. 
Does it hold to higher order? Empirically, factorization fails for the decays of P- 
waves into light hadrons at order a^. Why? 

2. Is there a rigorous field-theoretic definition of the long-distance factors? How 
do you compute |i?5(0)p and |i?p(0)p using lattice simulations? The long-distance 
factor |i?5(0)p was assumed to be the same for the decays of the ^5*0 and ^5"! states. 
It was also assumed to be the same for decays into light hadrons and decays into 
electromagnetic final states. At what point do these assumptions break down? 

3. How do you calculate relativistic corrections? Since ~ 1/3 for charmonium, 
the relativistic corrections may be more important than the next-to-leading order 
perturbative correction. For S-waves, one would expect the first relativistic correc- 
tion to have a long-distance factor that involves V^i?5(0), but this quantity diverges 
if the potential is Coulombic at short distances as in the case of QCD. How does 
one deal with this ultraviolet divergence? 

4. How do you take into account higher Fock states? While heavy quarkonium is 
predominantly in the Fock state \QQ), the amplitudes for higher Fock states, such 
as \QQg) and \QQgg), are certainly nonzero. At some level, they must contribute 
to annihilation rates. How can they be taken into account? 

These problems, which blocked progress on the theory of heavy quarkonium 
annihilation, resisted solution for more than a decade. Surprisingly, the solution 
to these problems is rather simple when viewed from the proper perspective. This 
perspective is provided by NRQCD, an effective field theory that makes it relatively 
easy to separate out effects of different momentum scales in heavy quarkonium. 
Once this framework is developed, the solutions to all the problems listed above fall 
into place. 

3. NRQCD and the General Factorization Formula 

3.1. Momentum Scales 

Heavy quarkonium is a system that involves several important momentum scales. 
First, there is the heavy quark mass M, which sets the length scale for the QQ 
annihilation process. Then there is the typical momentum Mv of the heavy quark, 
which sets the length scale for quarkonium structure. Next, there is Mv^, the typical 
kinetic energy of the heavy quark. This scale is important, because the effects of 
dynamical gluons with momenta of order Mv"^ can not be taken into account via 
an effective potential between the heavy quark and antiquark. They can however 
be taken into account through higher Fock states, such as \QQg). Finally, there is 
^QCDi the scale of nonperturbative effects associated with light quarks and gluons. 
It is the existence of so many important length scales that makes heavy quarkonium 



such an interesting and challenging problem. 

If the heavy quark mass M is large enough, than these momentum scales will 
be well-separated: 

Aqcd < {Mv^f < {Mvf < M^. (8) 

We will refer to 1/M as the short-distance scale, and the other 3 scales will be 
referred to collectively as long-distance scales. Our analysis of heavy quarkonium 
will be based on separating the effects of the short-distance scale from those of the 
long-distance scales. This separation is very distinct in bottomonium, where the 
typical value of v"^ is ~ 1/10, and it is reasonably clear in charmonium, where 
^2 ^ 1/3. 

3.2. NRQCD Lagrangian 

A convenient tool for separating the short-distance scale M from the long- 
distance scales in heavy quarkonium is nonrelativistic QCD (NRQCIDl. This is an 
effective field theory which was developed by Lepage and coUaboratorsSlll precisely 
for this purpose. It is designed to remove the scale M from the problem, while pre- 
cisely reproducing the bound-state dynamics at length scales of order l/(Mt>) or 
longer. The lagrangian for NRQCD is 

-^NRQCD = -Clight + -^^heavy + SC, (9) 

where Cnght is the usual relativistic lagrangian for the gluons and the light quarks. 
The heavy quarks and antiquarks are described by the term 

^heavy = (^tDo + ^ + (^^0 - ^) X, (10) 

where ip is a 2-component spinor field that annihilates a heavy quark, x is a 2- 
component spinor field that creates a heavy antiquark, and Z)^ is the gauge-covariant 
derivative. The lagrangian £iight + -Cheavy describes ordinary QCD coupled to a 
Schrodinger field theory for the heavy quarks and antiquarks. This field theory is 
nonrenormalizable, and requires an ultraviolet cutoff A whose natural scale is of 
order M. The final term 6C in the lagrangian (|^) includes all possible local gauge- 
invariant counterterms. Their coefficients must be tuned so that the ultraviolet 
divergences in NRQCD calculations of long-distance quantities are cancelled, and 
so that the residual finite terms reproduce the corresponding results from full QCD. 

The solutions to the heavy-quark field equation (z-Dq + D^/(2M))^/' = contain 
a dynamically-generated length scale Mv. This scale is generated by the non- 
perturbative balance between the effects of the operators Dq and D^/2M, with 
Dq ~ Mf^ and D ~ Mv. Thus there is also a dynamically-generated small pa- 
rameter V, which can be used to analyze the importance of the various terms in 
the NRQCD lagrangianll^. By using the field equations, one can show that the 
lagrangian Cnght + -^heavy can be used to calculate matrix elements between quarko- 
nium states up to an error of order v^. If higher accuracy is required, one must add 



the counterterms 
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The coefficients in (|TT]) can be determined by demanding that perturbative scat- 
tering amplitudes in NRQCD agree with those in full QCD to order v"^. Using the 
lagrangian £iight + -^heavy + ^^C., matrix elements between heavy quarkonium states 
can be calculated with errors of order w^. 

3.3. Annihilation in NRQCD 

The annihilation of a QQ pair cannot be described directly in NRQCD, be- 
cause annihilation processes such as QQ — > g* produce gluons with small wave- 
lengths of order 1/M. However, annihilation can be taken into account indirectly 
through its effects on QQ scattering amplitudes. At long distances, the amplitude 
for QQ — > — > QQ can be reproduced by a point-like interaction for QQ scatter- 
ing. The effects of annihilation can therefore be taken into account adding 4-fermion 
operators to the lagrangian: 

*A-,™.„^E^S#?^0„, (12) 

n 

where the sum is over all possible local 4-fermion operators On that annihilate and 
create a heavy-quark-antiquark pair, and dn is the scaling dimension of On- An 
example of such an operator is V'^xxV; which has a scaling dimension of 6. The 
coupling constants /„(as) can be computed by matching perturbative amplitudes 
for QQ scattering in NRQCD with the corresponding amplitudes in full QCD. 

According to the optical theorem, the annihilation rate of a QQ pair is pro- 
portional to the imaginary part of its forward scattering amplitude. In order to 
reproduce these imaginary parts in NRQCD, the coefficients /„ of the 4-fermion 
operators must have imaginary parts. These imaginary parts must also reproduce 
the effects of annihilation of heavy quarkonium. Up to a factor of (—2), the decay 
rate F can be identified with the imaginary part of the energy of the state. Since 
that energy is the expectation value of the interaction lagrangian, we can write the 
annihilation part of the decay rate of a heavy quarkonium state H as 

T{H) = 2 1m(ff|£4-fermion|^). (13) 

Using the expansion (p!2D, this becomes 



^ E ""ttf"V ia.|g>. (14) 

This is the general factorization formula for the annihilation decay rates of heavy 
quarkonium states0. It expresses the decay rate as the sum of long-distance matrix 



elements {H\On\H), multiplied by short-distance coefficients. The sum is over all 
possible 4-fermion operators. All the dynamics associated with the scale M resides 
in the short-distance coefficients, which are proportional to the imaginary parts of 
coupling constants in the NRQCD lagrangian. All effects from distance scales of 
order l/{Mv) or larger are contained in the matrix elements. 

3.4- Simplifying the Matrix Elements 

As it stands, the general factorization formula is not very useful because it 
contains infinitely many terms. To make it into a practical tool, we must exploit 
the dynamically-generated small parameter v in NRQCD. As f — 0, each of the 
matrix elements in (|I4D scales with a definite power of f . At any given order in 
V, only a finite number of the matrix elements contribute to the decay rate. The 
number of independent matrix elements is even smaller, because there are many 
relations between the matrix elements that hold at low orders of v. Below, we 
summarize the methods thatcan be used to reduce the number of terms in the 
general factorization formula.lii 

At a given order in f , the number of matrix elements can be reduced to a 
finite number by using velocity scaling rules for the matrix elements. These scaling 
rules consist of scaling rules for the operators and scaling rules for the probabilities 
of the Fock states that give the leading contributions to the matrix elements. The 
scaling rules for the operators are fairly simple Jl3 and can be summarized by D ~ f , 
Dq ~ v"^, g'Ei ~ f ^, and (^B ~ v^. These rules give the scaling behavior of the matrix 
element {H\On\H) if the operator On annihilates and creates the QQ pair in the 
dominant Fock state \QQ) of the meson \H). An example is the matrix element 
(^c|V'^XX^V'l'7c)5 where the operator creates and annihilates a cc pair in a color- 
singlet ^5*0 state. The matrix element {H\On\H) is suppressed by one extra power 
of if On annihilates and creates the QQ pair in the \QQg) Fock state that can 
be reached by an El transition from the dominant \QQ) state. An example is the 
matrix element {hc\ip'^T"'xX^T""ijj\hc) , where the operator creates and annihilates a 
cc pair in a color-octet ^5*0 state. In the dominant Fock state \QQ) of the he, the 
QQ pair is in a color-singlet ^Pi state, but the he also has an amplitude of order v 
to be in a \QQg) Fock state, with the QQ pair in a color-singlet ^5*0 state. Thus 
the matrix element is suppressed by a factor of f ^ from the probability of this Fock 
state. 

The matrix elements that contribute to a decay rate to a given order in v are 
not all independent. If the operator C„ creates and annihilates the QQ pair in the 
dominant Fock state \QQ), then the matrix element {H\On\H) can be simplified by 
using the vacuum-saturation approximation. For example, 

{r}cWxxH\Vc) = (^c|V'^x|0)(0|xVl^?c)(^l + 0{v^)y (15) 

The vacuum-saturation approximation provides relations between the matrix el- 
ements that contribute to decays into light hadrons and those that contribute to 



decays into electromagnetic final states. In NRQCD, this is a controlled approxima- 
tion with a relative error that scales like v^. There are also relations between matrix 
elements that follow from heavy-quark spin symmetry. This is an exact symmetry 
of the lagrangian £iight + -Cheavy Thus the relations implied by heavy-quark spin 
symmetry are correct up to errors of relative order f^. An example of a relation 
that follows from heavy-quark spin symmetry is 

(OlxVl^c) = e-{Q\xH\i'{e))i^ + 0{v^)y (16) 

where e is the polarization vector of the if). 



4. Applications 

S-waves at Leading Order in v"^ 
In previous calculations for S-wave states, the annihilation decay rates were as- 
sumed to factor into short-distance coefficients multiplied by |-Rs(0)p, where Rs{0) 
is the radial wavefunction at the origin. This long-distance factor was assumed to 
be the same for decays of the 0"+ state and the 1 state. It was also assumed to 
be the same for decay rates into electromagnetic final states and for inclusive decay 
rates into hadrons. Our general factorization theorem verifies that these assump- 
tions are correct at leading order in f^. At leading order in f^, the decay rates of 
the ^ state rjc and the 1 state rjc can be written 
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In (|T^) and (|T§), the vacuum-saturation approximation has been used to replace the 
matrix elements (^d'^^XxVl^c) and {ip\ip'^ crx-X^ ^i^li^) by the same factor \Rs\'^ that 



appears in the electromagnetic decay rates (^ and (|20|) . The radial wavefunctions 
at the origin for the ip and rjc are given by the matrix elements 



Rr, 




R^ = ^e-(0|xWlV^(e))- 



(21) 
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At leading order in f^, we can use heavy-quark spin symmetry to equate Rs in (p!7|)- 
([20| ) with Rri^, or with i?^, or with their weighted average: Rs = {R^^ + 3-R^)/4. 



At this level of accuracy, the only advantage of the NRQCD formalism is that it 
provides a rigorous nonperturbative definition of the long-distance factor l-^ii, so 
that it can be calculated with lattice simulations. It has been calculatedtj'Ej to 
leading order in v"^. The short-distance coefficients Fi in the factorization formu- 
las (P!7D-(|2I1|) are proportional to the imaginary parts of coupling constants in the 
NRQCD lagrangian. They have all been calculated to next-to-leading order in a^. 

4-2. S-waves at Next-to-leading Order in v"^ 

At next-to-leading order in f^, the factorization formula for S-wave states be- 
comes less trivial, with 3 independent long-distance factors: 
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In the first terms on the right sides of (p3|)-(|26|), the radial wavefunctions at the 
origin i?^^ and R^ can not be equated with their weighted average Rs without loss 
of accuracy. The third independent matrix element V^Rs can be interpreted as a 
renormalized laplacian of the wavefunction at the origin: 
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The short-distance coefficients in (|23| )-(P6D are the same as in (|T7D-(|20|), and 
are known to next-to-leading order in a^. The short-distance coefficients Gican be 
extracted at leading order in as from calculations by Keung and Muzinichll^, who 
developed a phenomenological model for taking into account relativistic corrections 
to quarkonium decays. A phenomenological analysis of S-wave decays, including 
the relativistic corrections, is in progresstZl. 

4-3. P-waves at Leading Order in v"^ 

The correct factorization formulas for decays of P-wave states! are nontrivial 
even at leading order in f^: 
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There are two independent long-distance factors that appear in these decay rates. 
The long-distance factor |-Rpp can be identified with the square of the derivative of 
the wavefunction at the origin for the or any of the other P-states: 



R'p ^ ^e-{0\x^V^\Kie)) 



(31) 



Coulomb 



The second long-distance factor in ([ 
operator in NRQCD: 



and (l29|) is the matrix element of a 4-fermion 

(32) 



The operator in the matrix element annihilates and creates a QQ pair in a color- 
octet S-wave state, and thus measures the probability density for the QQ pair 
to be in a color-octet S-wave state at the origin. This term in the factorization 
formula represents a contribution to the decay rate from the \QQg) component of 
the quarkonium wavefunction. The matrix element (Og) depends logarithmically 
on the ultraviolet cutoff A of NRQCD, which can be interpreted as an arbitrary 
factorization scale in the full theory: 
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The A-dependence of (Os) must be cancelled by the coefficients Fi in the factoriza- 
tion formulas (|28|) and (^9]). The logarithmic dependence of these coefficients on A 
explains the infrared divergences discovered by Barbieri et al. in their pioneering 
calculations of P-wave decay sii. 

The short-distance coefficients Fi(3Po), i^i(^P2), ^8(^5'o), and F^{^Si) in (||)- 
(|30| ) are knownl to leading order in while F^^(^Po) and F^^(^P2) have been cal- 
culated to next-to-leading orderB. Accurate lattice calculations of the long-distance 
factors are not yet available. Fortunately the long-distance factors can be deter- 
mined phenomenologically. Experiment E760 at Fermilab^, which studies the reso- 
nant production of charmonium in collisions, has provided precise measurements 
of the total widths of the Xci and Xc2"- They can be used to determine R'p and (Cs), 
and then the remaining four P-wave decay rates can be predicted. The predictiiin 
for Xc2 ~* 77 is consistent with recent measurements by E76Clli3 and by CLE0c2l. 
The predictions for the total widths of the he and Xco and for the rate for Xco 77 
should be tested in the foUowup experiment E835. 



5. Summary 

A rigorous theory for the annihilation decay rates of heavy quarkonium has 
recently been developed. The inclusive decay rate of a quarkonium state H into 
light hadrons satisfies a general factorization formula: 



(34) 



The sum is over all 4-fermion operators in NRQCD that create and annihilate 
both a quark and an antiquark, but at any given order in f^, there are only a 
finite number of operators that contribute. The long-distance factors {H\On\H) 
are well-defined matrix elements that can be calculated numerically using lattice 
simulations of NRQCD. The short-distance coefficients can be calculated as 
perturbation expansions in as{M). Using this general factorization theorem, it is 
now possible to calculate the annihilation decay rates of heavy quarkonium from 
first principles, with the only inputs being the heavy-quark mass M and the QCD 
coupling constant. 
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